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ABSTRACT 

An unbiased method for improving the resolution of astronomical images is presented. 
The strategy at the core of this method is to establish a linear transformation between 
the recorded image and an improved image at some desirable resolution. In order to 
establish this transformation only the actual point spread function and a desired point 
spread function need be known. Any image actually recorded is not used in establishing 
the linear transformation between the recorded and improved image. 

This method has a number of advantages over other methods currently in use. It 
is not iterative which means it is not necessary to impose any criteria, objective or 
otherwise, to stop the iterations. The method does not require an artificial separation 
of the image into "smooth" and "point-like" components, and thus is unbiased with 
respect to the character of structures present in the image. The method produces 
a linear transformation between the recorded image and the deconvolved image and 
therefore the propagation of pixel-by-pixel flux error estimates into the deconvolved 
image is trivial. It is explicitly constrained to preserve photometry and should be 
robust against random errors. 

Key words: methods : data analysis - methods : numerical - techniques : image 
processing 



1 INTRODUCTION 

In astronomy the problem of correcting images for imper- 
fect telescope optics, atmospheric turbulence and other ef- 
fects that adversely influence the image quality is very well 
known. There exist therefore many different ways to improve 
on images or imaging techniques in order to obtain more 
detailed spatial information. First of all there are hardware 
based solutions such as adaptive optics and interferometry. 
It is of course always desirable to make use of such tech- 
niques whenever possible. However, given a recorded image 
with a known point spread function (PSF) and an estimate 
of the point by point error it is still possible to improve 
on the spatial resolution by means of software : numerical 
image deconvolution. 

A number of algorithms already exist that attempt 
to do this image reconstruction. Best known in the field 
of interferometry at radio wavelengths are probably the 
CLEAN method (Hogbom, 1974 ; Schwarz, 1978 ; Wakker & 
Schwarz, 1988) and the maximum entropy method (MEM) 
(cf. Narayan & Nityananda, 1986). For the deconvolution 
of optical images the Richardson-Lucy (RL) method is well 
known (cf. Richardson, 1972 ; Lucy, 1974, 1992, 1994). Quite 
recently a new method was presented by Magain, Courbin, 
& Sohy (MCS, 1998). A characteristic of all of these meth- 
ods is that images are reconstructed by placing flux (MEM 
and RL), or building blocks such as point sources (MCS, 



CLEAN, two channel RL) in the field of view and minimiz- 
ing the difference between the image of this model (after 
convolution with the PSF) and the actual image. Since this 
is an inverse problem the solution is generally not unique 
and it can be quite sensitive to errors in the data. Thus in 
the minimization there is some need for regularization which 
is usually a smoothness constraint. The various method dif- 
fer in which building blocks are used and in the form of 
regularization applied. 

It is however somewhat unsatisfactory to proceed in this 
manner since astronomical images are not generally easily 
described by just point sources or objects of a certain shape 
or size, and may not conform to the smoothness constraint 
applied. Forcing an algorithm to nevertheless build the im- 
age up within such constraints may well introduce an un- 
desirable bias. For this reason it is useful to consider an 
alternative that does not assume anything about properties 
of the image in the deconvolution. In fact the method pre- 
sented here does not even use the image actually recorded 
until its very last step. 

The method of subtractive optimally localized averages 
(SOLA) was originally developed for application in the field 
of helioseismology (cf. Pijpers & Thompson, 1992, 1994) to 
determine internal solar structure and rotation. Since then 
it has also been successfully applied to the reverberation 
mapping of the broad line region of active galactic nuclei 
(AGN) (Pijpers & Wanders, 1994). Instead of operating on 
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the image, the SOLA method uses the PSF with which the 
image was recorded. With this PSF and a user-suppHed de- 
sired PSF a Unear transformation is constructed between 
any recorded image for which that PSF apphes and its de- 
convolved counterpart. The resolution that can be attained 
in this way is only limited by the sampling of the recording 
device (the pixel size of the CCD) and by the level of the 
flux errors in the recorded image. Since the transformation 
is linear, it is quite straightforward to impose photometric 
accuracy. Astrometric accuracy at the pixel scale is simi- 
larly guaranteed since there is no 'positioning of sources' in 
the image by the algorithm. Sub-pixel accuracy, claimed for 
some deconvolution methods, implies subdividing each pixel 
into subpixels. It requires knowledge of the PSF at very high 
accuracy and very small errors in the data in order to de- 
convolve down to a sub-pixel scale. If such information is 
available the SOLA method can easily accommodate sub- 
pixel scale deconvolution, without substantial modifications. 
In what follows however, it is assumed that a single pixel is 
the smallest scale required. 

In section 2 the SOLA method is presented. In section 3 
the method is applied to an example image to demonstrate 
the workings of SOLA. In section 4 the method is applied to 
some astronomical images. Some conclusions are presented 
in section 5. 



2 THE SOLA METHOD 
2.1 arbitrary PSFs 

The strategy of the SOLA method in general is to find a set 
of linear coefficients c which, when combined with the data, 
produce a weighted average of the unknown convolved func- 
tion under the integral sign, where the weighting function is 
sharply peaked. In the application at hand this means find- 
ing the linear transformation between an image recorded at 
a given resolution and an image appropriate to a different 
(better) resolution. 

The relation between a recorded image D and the actual 
distribution of flux over the field of view / is : 

D{x,y) = j dx'dy' K{x',y' ;x,y)I{x',y') (1) 

where K is the PSF. If one assumes that the PSF is constant 
over the field of view then : 



K{x', y' ■ X, y) = K{x - x' ,y - y') 



(2) 



Of course generally D is not known as a continuous function 
of {x,y), but instead it is sampled discretely as for instance 
an image recorded on a CCD. Thus one has as available data 
the recorded pixel- by-pixel values of flux D(xi,yj). These 
measured fluxes will usually be corrupted by noise and thus 
the discretized version of equation (1) is : 



(3) 



Dij = j dx'dy' Kij{x',y')I{x',y') +nij 

where now Kij refers to the PSF appropriate for the pixel 
at {xi,yj) and Dij is the flux value recorded in that pixel. 
In the vocabulary usual for the SOLA method the Kij are 
referred to as integration kernels. 

In the SOLA technique a set of linear coefficients c; 
is sought which, when combined with the data, produces a 



value for the flux R in any given pixel that would correspond 
to an image recorded with a much narrower PSF. Writing 
this out explicitly and using (3) yields : 



j dx'dy' ^^^ciKi{x',y ')}l{x',y') + Y^ 



(4) 



in which the double subscript ij has been replaced by a 
single one I for convenience. Thus one would construct the 
c; such that the averaging kernel K, deflned by : 



K,{x',y') = y^ciKi{x',y') 



(5) 



is as sharply peaked as possible. If one does this for all lo- 
cations on the CCD the collected values Rm are then the 
fluxes corresponding to the image at this (better) resolution 
with a (improved) "point spread function" /C. The so-called 
propagated error, the error in the flux R is : 



(6) 



Here the 7V;„i is the error variancc-covariance matrix of the 
recorded CCD images where both / and m run over all {i, j) 
combinations of the pixel coordinates. If the errors are un- 
correlated between pixels then (6) reduces to : 

o-fl^^Qo-f (7) 

which is trivially computed once the coefficients c; are 
known. 

Ideally one would wish to construct an image corre- 
sponding to an infinitely narrow PSF : a Dirac delta func- 
tion. In practice this cannot be achieved with a finite amount 
of recorded data. As has already been pointed out by Magain 
et al.(1998) one must in the deconvolved image still satisfy 
the sampling theorem. A further restriction arises because 
of the noise term in equation (3) . As is well known in helio- 
seismology the linear combination of data corresponding to 
a very highly resolved measurement usually bears with it a 
very large propagated error. In order to obtain a fiux value 
for each pixel in the deconvolved image that does not have 
an excessively large error estimate associated with it, one 
needs to remain modest in the resolution sought for in the 
deconvolved image. 

Finding the optimal set of coefficients taking these limi- 
tations into account can be expressed mathematically in the 
following minimization problem. One needs to minimize for 
the coefficients c; the following : 



J dxdy [K, - Tf + ji CiC„ 



.Ni, 



(8) 



Here /i is a free parameter which is used to adjust the relative 
weight given to minimizing the errors in the deconvolved 
image and to producing a more sharply peaked kernel K,. 
The higher the value of fi the lower this error but the less 
successful one will be in producing a narrow PSF. In SOLA 
one is free to choose the function T. A common choice in 
SOLA applications is a Gaussian : 



T = 



■ exp 



{x - xof + {y- yof 
A2 



(9) 



Here (xo,2/o) is the location for which one wishes to know 
the fiux at the resolution corresponding to the width A. / 
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is a normalization factor chosen such that 
dxdy T=l 



I 



(10) 



although any set of locations {xo,yo) can be chosen, a nat- 
ural choice in the application at hand is to take all original 
pixel locations {xi,yj). If one wishes to deconvolve to sub- 
pixel scales this can be done by an appropriate choice of the 
{xo,yo) and A. 

In terms of an algorithm the problem of minimizing the 
function (8) leads to a set of linear equations : 



The elements of the matrix A are given by 
Aim = J dxdy Ki{x,y)Km{x,y) + i^Nim 
The elements of the vector b are given by : 
dxdy T{x,y)Km{x,y) 



(11) 



(12) 



(13) 



Writing out the dependencies on the free parameters explic- 
itly, determining the coefficients ci results from a straight- 
forward matrix inversion : 



ci{xo,yo ■,A,ii) = AiJ^{fx)bm{xo,yo ;A) 



(14) 



It is clear that for each point (xo,J/o) there is a separate set 
of coefficients q which will depend on the resolution width 
A required and on the error weighting /x. Note that it is not 
necessary to invert a matrix for every location (xo, yo), which 
would certainly be prohibitive if one wishes to calculate the 
entire deconvolved image. For a given error weighting fx one 
needs to invert A only once. Only the elements of the vec- 
tor b need be recomputed for different locations or different 
resolutions. 

In order to ensure that at every point in the recon- 
structed image the summed weight of all measurements is 
equal and thus a true (weighted) average it is necessary to 
additionally impose the condition : 



(15) 



It is this condition that imposes photometric accuracy on the 
reconstructed image. Using the method of Lagrange multi- 
pliers this condition is easily incorporated into the matrix 
equation (11) by augmenting the matrix A with a row and 
column of I's, and a corner element equal to 0. The vector 
b gains one extra element equal to 1 £is well. The details of 
this procedure can be found in Pijpers & Thompson (1992, 
1994). 

2.2 translationally invariant PSFs 

Although the method described above can work in princi- 
ple with general PSFs K, the matrix inversion becomes in- 
tractable very quickly as the number of pixels increases. For 
an image of M x M pixels the number of elements in the 
matrix A is x M^. The matrix A is symmetric but even 
so a naive matrix inversion routine would require a number 
of operations scaling as M^. 

However the entire procedure for obtaining the trans- 
formation coefficients for all locations of the CCD can be 



speeded up considerably if one accepts some restrictions for 
the properties of the PSF K and of the expected errors Nim- 
The first restriction is to assume that the PSF is constant 
over the field of view, that is to say that equation (2) is 
valid. When condition (2) is met one can easily demonstrate 
that in equation (12) the integrals of the cross products of 
the PSFs are a convolution : 



/ 



dxdy Ki{x,y)Km{x,y) 



I 



dxdy K{x -Xi^,y- yj^ )K(x -Xi^,y- yj^ ) (16) 



= Jdx'dy'Kix',y')K'iAx.^.,-x' 



in which : 
K'{x -x ,y- 
Axi 



K(x' ■ 



x,y 



x',Ayj^ji -y 



V) 



(17) 



h = Vim. - Vh 



Evaluating all the elements of the matrix A is much 
simplified by doing this two-dimensional convolution as a 
multiplication in the Fourier domain. This calculation is 
then dominated by the FFT calculation which requires 
C'(M^log(M)) operations. Similarly the vectors 6 in (13) 
can be evaluated for all locations {xi,yj) with a single two- 
dimensional convolution of K and T, again dominated by 
the FFT. 

If the CCD pixels are assumed to be equally spaced the 
matrix A for /j, = can be constructed in such a way that it 
becomes of a special type known as symmetric block circu- 
lant with circulant blocks (BCCB) , for which very fast inver- 
sion algorithms exist. Circulant matrices have the property 
that every row is identical to the previous row, but shifted to 
the right by one element. The shifting is 'wrapped around' 
so that the first element on each row is equal to the last ele- 
ment of the previous row. Thus the main diagonal elements 
are all equal and on every diagonal parallel to the main di- 
agonal of the matrix all elements are equal as well. A BCCB 
matrix is a matrix that can be partitioned into blocks in such 
a way that each row of blocks is repeated by shifting (and 
wrapping around) by one block in the subsequent row of 
blocks and each individual block is circulant. It can be shown 
that circulant matrices can be multiplied and inverted using 
Fourier transforms, and by extension BCCB matrices can 
be multiplied and inverted using two-dimensional Fourier 
transforms. The detailed steps of the algorithm are worked 
out in the appendix. 

The restriction on the matrix Nim is that is must also 
be a symmetric BCCB matrix for the fast inversion algo- 
rithm to work. It is evident that fully optimal results can 
only be obtained if the full N'^ x JV^ covariance matrix of 
the errors is used. However, the error correlation function for 
the pixels is expected to behave similarly to the point spread 
function in the sense that it is large (in absolute value) for 
small pixel separations and small for large pixel separations, 
independently of where on the CCD the pixel is located. It is 
therefore likely that the error covariance matrix will already 
be BCCB or be very nearly so. Since its role in the minimiza- 
tion of (8) is to regularize the inversion it is in practice not 
essential that the exact variance-covariance matrix be used. 
Experience in using SOLA in other fields has shown that the 
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Figure 1. The 128 X 128 pixels image used in testing the algorithm. Top left panel ; the original image. Top right panel : the original 
convolved with the target PSF ; a Gaussian with A = 1.5 pixels. Bottom left panel : the original convolved with a PSF which is the 
sum of a Gaussian with A = 10 pixels and an 0.1% contribution from a Gaussian with A = 1 pixel. Bottom right panel : the image 
after SOLA deconvolution of the bottom left image. In all images the grey-scale is linear. In the bottom left image noise is added before 
deconvolution. In the bottom right image the noise propagated in the deconvolution has an expectation value of ~ 0.5 in arbitrary flux 
units and the S/N ratio for the brightest pixel is ~ 1000. 



results of linear inversions are robust to inaccuracies in the 
error matrix, as long as those are not orders of magnitude 

large : if for instance substantial amounts of data (fluxes in 
pixels) are to be given small weight in the resulting linear 
combination, because of large errors associated with them, 
this can give rise to large departures from BCCB behaviour 
of the error covariance matrix. This would then cause prob- 
lems for the fast version of the SOLA method presented here. 
Thus if Nim is not circulant it should in most cases be sufH- 
cient to use a BCCB matrix that is close to the original : one 



could think of using a modified matrix A'^;^ on the diagonals 
of which are the average values over those diagonals of the 
true Nim. Of course once the coefficients have been deter- 
mined, when calculating the propagated errors one should 
use equation (6) with the proper variance-covariance matrix 

Nlm. 

As is shown in the appendix the matrix corresponding 
to the collection of all vectors of coefficients c, which results 
from the multiplication of A with the matrix corresponding 
to the collection of all vectors 6 (one for every pixel) , is also 
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a BCCB matrix. The process of combining these coefficients 
Q with the recorded fluxes on the CCD to form the improved 
image is : 



Ri 



CkiDi+k-ij+i-i 



(18) 



k I 



Since the matrix C is a BCCB matrix and therefore its trans- 
pose is as well, the following holds : 

Rij = y^ CklDj+k-lj+l-l = y^ CjkPj+k-l j+Z-l 



k, I 



2-k2-lDi+k-l j+l-l 



(19) 



Ck' I'Di+i-k' j+i- 



k', I 



Prom the final equality in (19) it is clear that the process 
of combining the matrix of coefficients with the image is a 
convolution, and hence can also be done using FFTs. 

From the above it is clear that limiting the algorithm 
to the case of a PSF that is constant over the CCD implies 
a profound reduction of the computing time. If the PSFs K 
satisfy the condition (2) the vectors b collected together for 
all {xo,yo) ~ {xi,yj) form a BCCB matrix, and therefore 
the matrix inversion of A and its subsequent multiplica- 
tion with all vectors b, shown in equation (14), can be done 
in 0{M^ log(M)) operations. The entire deconvolved image 
is thus produced in 0{A4'^ log(M)) operations. This accel- 
eration of the algorithm over the version described in the 
previous section is so substantial that even when the PSF 
is not constant over the CCD it is worthwhile subdividing 
the image into subsections in which the PSF can be closely 
approximated by a single function K. The error introduced 
in this way can be estimated in a way similar to what is done 
in the application of SOLA to the reverberation mapping of 
AGN (Pijpers & Wanders, 1994), and should generally be 
much smaller than the propagated error from equation (6). 
If such a subdivision is undesirable, there is the possibility 
of reverting to the more general algorithm of section 2.1. 
For a single peaked PSF the matrix A should have a banded 
structure to which fast sparse matrix solvers can be applied. 
In this case one could use the inverse of the matrix A for 
an approximated PSF that is translatioually invariant as a 
pre-conditioner to speed up the matrix inversion for the case 
of the true PSF. 



3 APPLICATION TO A TEST IMAGE 
3.1 constructing a narrow PSF 

In the first instance it is useful to test the algorithm on a test 
image for which the result and the errors are known. To this 
end an artificial image of a cluster of stars is convolved with 
two different PSFs. One PSF is a sum of two Gaussians ; 
one with a width A = 10 pixels in which 99.9% of the flux is 
collected, and a second one with a width A = 1 pixels which 
collects the other 0.1% of the flux. Poisson distributed noise 
is added to every pixel and this 'dirty' image serves as the 
image to be deconvolved. The other PSF is a Gaussian with 
a width A = 1.5 pixels which is also the target chosen for 
the SOLA algorithm. Thus the deconvolved image can be 
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Figure 2. The difference between the image that is SOLA decon- 
volved and the original image convolved to the target PSF. The 
gray scale is adjusted so that the full scale is 0.01 X the scale in 
the right-hand side images of figure 1, which corresponds to lOcr 
of the noise in the deconvolved image. 
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Figure 3. A slice through the peak of the target PSF specified 
in the SOLA algorithm (solid line) and the averaging kernel K 
(dashed) constructed from the linear combination of the pixel 
PSFs. 



compared directly with the image obtained from direct con- 
volution of the original with the narrow target PSF. The 
results are shown in figure 1. In order to get an optimal re- 
production of the target PSF, the error weighting parameter 
/u is chosen to be equal to 0. 

It is clear that the bottom and top right panels are very 
similar and thus the image appears to be recovered quite 
well. To illustrate this further the two images can be sub- 
tracted. Figure 2 shows the SOLA deconvolved image minus 
the image convolved with the target PSF, with an adjusted 
gray scale to bring out the differences, which in the cen- 
tral portion of the image are all < 1%. Although there is 
no strong evidence for it in this image, the deconvolution 
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Figure 4. The error magnification A as a function of the width 
A of the target PSF specified. The PSF of the blurred image is 
as described in the text in all cases, the error weighting is /* = 0. 

can suffer from edge effects because part of the original im- 
age can 'leak away' in the convolution with the broad PSF. 
When deconvolving, the region outside the image is assumed 
to be empty and so a spurious negative signature is then in- 
troduced in the image. The magnitude of such edge effects 
must clearly depend both on the image and on the PSF of 
the 'dirty' image, since they are determined by the informa- 
tion that has been lost at the edges of the CCD. Of course it 
is desirable to demonstrate this method on a more realistic 
suite of images than just the simple one used here, which is 
work currently in progress. 

The averaging kernel that is constructed cannot in gen- 
eral match perfectly the target form, even in the absence 
of errors. In general any function can be completely recon- 
structed only out of a complete set of base functions. Since 
function space is infinite dimensional this would require an 
infinite number of base functions. In this test image there 
are no more available than the 128^ PSFs corresponding to 
each of the pixels and so there can never be a perfect match- 
ing of K. with T. In figure 3 a section through the maximum 
of both T and IC is shown. It is clear that at the 10 ppm 
level, the constructed averaging kernel starts getting wider 
than the target. If the ratio of the widths of the target form 
and actual PSF is even smaller than for this image, alter- 
nating negative and positive side lobes can show up in the 
averaging kernel which cause ringing. The amplitude of the 
sidelobes, and the width A below which ringing starts oc- 
curring, will in general depend on the weighting of the errors 
/Lt, as has been demonstrated from the application of SOLA 
to helioseismology (Pijpers & Thompson, 1994). 

As it stands the SOLA algorithm does not impose pos- 
itivity on the image. One could attempt to use a positiv- 
ity constraint to extrapolate the image beyond the recorded 
edges in such a way that it eradicates any negative fluxes in 
the image, which would in principle also remove associated 
positive artifacts around the edges. However, in the presence 
of errors this might be somewhat hazardous. Furthermore, in 
the presence of errors any edge effects might well disappear 
into the noise. 

If one assumes that the covariance of errors between 
pixels is equal to and the flux error in each pixel is equal 
to CT^ , or A'' = cr^J, then it is particularly simple to calculate 
the flux error for each pixel in the deconvolved image from 




Figure 5. The difference on an arbitrary magnitude scale be- 
tween the magnitude of the stars in the deconvolved image mdocon 
and the magnitude maim of their counterparts in the reference im- 
age constructed by convolving the original with the target PSF. 
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Figure 6. The position difference in pixel units between the 
stars in the deconvolved image mdecon and their counterparts 
in the image constructed by convolving the original with the tar- 
get PSF. Open circles : all stars with magnitudes between 15.9 
and 16.9, crosses : stars with magnitudes between 16.9 and 17.9, 
open squares : stars with magnitudes between 17.9 and 18.9, open 
triangles ; all stars with magnitude greater than 18.9. 



equation (7) since it is 
2 \ ^ 2 _ 



2 



A 2 2 
A (7 



(20) 



This factor A is usually referred to as the error magnifica- 
tion and is equal for all pixels in the reconstructed image. In 
general A increases as the ratio of Atarget / Apsf decreases. 
For the deconvolved image of figure 1 the error magnifica- 
tion is ^ 321. The magnitude of the error magnification 
for this simple example illustrates that the true limitation 
of deconvolving images may in practice not lie in the sam- 
pling theorem, but instead in the S/N of the recorded image. 
For example for a point source the peak flux in its central 
pixel will increase in the deconvolution by a factor which is 
roughly FWHMpsF /FWHMtargct . The noise will increase by 
a factor A and so the signal-to-noise ratio for point sources 
will scale roughly as : 



1_ / FWHMpsF \ f S 

A \^ FWHMtargct 



(21) 



Thus in the example shown the signal-to-noise ratio for point 
sources is degraded by a factor of roughly ~ 48 between the 
dirty image and the deconvolved image. This clearly requires 
a very high a signal-to-noise ratio in the dirty image, which 
means that in practice as dramatic a resolution enhancement 



as attempted here will not usually be possible. A more mod- 
est resolution enhancement of around a factor of 2 should 
in most cases be possible however, as can be deduced from 
figure 4. 

In order to show the relation between resolution and 
error magnification in figure 4 is shown the value of the 
error magnification as a function of the width A specified for 
the target function. If only the broad component had been 
present the error magnification would have boon unity for a 
target A — 10. Efi'ectively because of the narrow component 
which captures a mere 0.1% of the flux the image can be 
deconvolved to a PSF with A = 8 pixels without significant 
penalty in the magnification of the errors. 

The CPU time used to construct the matrix A, all the 
vectors b, to invert A and multiply with all b and finally to 
combine the coefficients with the image takes ~ 0.5 min on 
an SGI workstation for this 128 x 128 image. 

3.2 photometric and astrometric accuracy 

Since there is no placement of flux or point sources, the al- 
gorithm should automatically be astrometrically accurate. 
Photometric accuracy is ensured by explicitly constrain- 
ing the linear coefficients to sum to unity, i.e. imposing 
constraint (15). In order to demonstrate both these prop- 
erties for this test image a standard photometric package 
DAOPHOT was used to do aperture photometry on the de- 
convolved image, and on the image obtained by convolving 
the original with the target PSF. The deconvolved image of 
course has noise propagated from the dirty image. The other 
image is kept noise-free to properly serve as a reference. The 
errors arc calculated by DAOPHOT and arc consistent with 
what is expected from the noise in the deconvolved image. 
In figure 5 is shown the difference between the magnitudes 
of the stars in the two images as a function of the stellar 
magnitudes in the reference image. The difference is clearly 
consistent with zero over the entire range of 5 magnitudes, 
and does not show any trend. 

DAOPHOT calculates the error bax assuming that the 
error in different pixels is uncorrelated. For the SOLA de- 
convolved image this is not the case. In the deconvolved 
image the error correlation function falls below 0.01 in ab- 
solute value only at inter-pixel distances larger than 12 
for this test case. Furthermore, because fluxes are combined 
with positive and negative coefficients, the error is not dis- 
tributed as for a Poisson process. If this is taken into account 
properly the error bars in figure 5 should be decreased and 
are then compatible with the actual scatter of the points. 

DAOPHOT also determines the positions of point 
sources in the image and therefore those positions can 
be used to determine astrometric accuracy. In figure 6 
are shown the difference in units of a pixel between the 
DAOPHOT determined positions in the two images. Here 
the stars have been grouped into 4 magnitude bins, each 
bin 1 magnitude in range, and starting from the brightest 
star with magnitude 15.9. The right-hand panel is a blow- 
up of the central portion of the left-hand panel. Figure 4 
shows that there is a trend in that the fainter stars show 
a greater scatter in position, the largest position difference 
being of the order of ~ 0.2 pixels. In the right-hand panel 
it can be seen that for stars brighter than magnitude 19 
the difference in positions is smaller than 0.03 pixels. These 
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Figure 7. Images of the galaxy UGC 5041. The top image is con- 
volved with the same broad PSF used on the image at the bottom 
left in figure 1. The middle image is the deconvolved image with 
a PSF with a FWHM of 2.5. The bottom image is the difference 
between the original image convolved with the target PSF, and 
the deconvolved image, no noise has been added. The gray scale 
of the bottom panel extends between ±0.1% of the gray scale of 
the middle image. 
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Figure 8. Image of a detail of the Orion Molecular Cloud obtained using adaptive optics in IR lines of shocked molecular hydrogen, 
north is at the top of the image, east is to the left. The image size is 256 X 256 pixels at 50mas/pixel. The top left image is as obtained 
with the Adonis instrument on ESO's 3.6iu telescope. The PSF for this image is shown on a Ix enlarged scale as the 32 X 32 pixel image 
in the top left-hand corner of the top right-hand panel. The bottom left- and right-hand panels are the deconvolved image using different 
gray-scales. The two bottom images in the top right-hand panel show 32 x 32 pixels images of the PSF for the deconvolved image, using 
the same gray scales as the corresponding images, and the same spatial scale as the PSF for the original. For all images the dynamic 
range between lightest and darkest colour is a factor of ~ 8 in flux level. 



uncertainties arc entirely consistent with the accuracy with 
which DAOPHOT can determine stellar positions. 

Prom figures 5 and 6 it is clear that if any errors in pho- 
tometry or in position arc introduced by the deconvolution 
process, they are much smaller than the errors due to the 
random noise. 



4 APPLICATION TO ASTRONOMICAL 

IMAGES 

4.1 UGC 5041 

To give a somewhat more interesting example, a high res- 
olution HST image of a galaxy is blurred and then decon- 
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volved to demonstrate that the method also works on an 
image which contains a combination of extended structure 
and point sources. 

The galaxy UGC 5041 is an Sc type galaxy at a rodshift 
of 0.027 (Hayncs et al.., 1997). It has been part of various 
•surveys for use in studies of clustering and in establishing 
distance scales for the Tully-Fisher distance method. In fig- 
ure 7 is shown a 512 x 512 image of this galaxy obtained in 
March 1997 using the WFPC2 (WF3) instrument on board 
the Hubble Space Telescope (HST) with the F814W filter 
which corresponds to the I-band. The resolution of the orig- 
inal image has a FWHM of 1.4 pixels and is convolved with 
the same broad PSF used in the bottom left panel of fig- 
ure 1. It is then deconvolved to the same resolution as the 
right-hand images of figure 1 and the result is shown in the 
middle panel of figure 7. The difference image between this 
deconvolved image and a reference is shown in the bottom 
panel of figure 7, where the gray scale is enhanced to demon- 
strate that in the absence of noise the differences between 
the deconvolved image and the reference image are less than 
0.1% of the peak flux apart from edge effects. 

4.2 Orion Molecular Cloud 

Observations in IR lines of shocked H2 of the SE part of the 
Orion Molecular Cloud complex (OMCl) have been per- 
formed at the ESO 3.6m telescope taking advantage of the 
high spatial resolution given by adaptive optics (Adonis at 
50 marcsec/pixel) combined with the high spectral resolu- 
tion given by a Fabry-Perot (R=1000) (Vannier et al., 1998). 
The image is deconvolved using a target PSF with FWHM 
3 pixels and an error weighting /x = 0. The resulting error 
magnification factor is A = 6.2. The dynamic range in the 
deconvolved image is ^ 30 as opposed to ~ 8 in the origi- 
nal. For this reason the deconvolved image is shown twice : 
bottom left in figure 8 is shown a linear gray scale extending 
from an estimated noise level to 8 times that, bottom right 
is shown a linear gray scale extending from the peak level 
in the deconvolved image which is ~ 30 times the estimated 
noise level, to 1/8 of that. The PSF of the deconvolved im- 
age, i.e. the constructed averaging kernel K., is shown using 
the same two gray scales in the bottom part of the top right- 
hand panel of figure 8. 

Although some of the 'graininess' in the bottom panels 
must be due to the increased noise compared to the top left- 
hand image, fine structure can clearly be seen in the bottom 
panels of figure 8. There is also some evidence of edge effects 
at the top and right of the image. From the image of the PSF 
of the deconvolved image it is also clear that, as expected, 
the Gaussian target is not reproduced perfectly over the 
entire dynamic range of the deconvolved image. Comparing 
the PSF images with the same dynamic range from the peak 
down (top left and bottom right in the top right-hand panel 
of figure 8) it is clear that the PSF is indeed much narrower 
for the deconvolved image. 

5 CONCLUSIONS 

In this paper the SOLA inversion method, well known in he- 

lioseismology, is applied to the reconstruction of astronomi- 
cal images. It is demonstrated how a linear transformation is 
constructed between any image recorded with a known PSF 



and its deconvolved counterpart with a different (narrower) 
PSF. The method itself uses only the PSF and no assump- 
tions are made concerning what is contained within the im- 
age(s) to be deconvolved. It is furthermore shown that in 
the case of translationally invariant PSFs, a fast algorithm, 
using 0{N log N) operations where A'^ is the total number 
of pixels in the image, can be constructed, which allows de- 
convolution of even 1024 x 1024 images within half an hour 
on medium-sized workstations. 
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APPENDIX A: THE INVERSION AND 
MULTIPLICATION OF CIRCULANT 

MATRICES 

As it stands the matrix described in (16) does not conform 
to the criteria for a block circulant with circulant blocks 
(BCCB) matrix. Instead it is a block Toeplitz matrix with 
Toeplitz blocks. For the latter type of matrix the elements 
are also identical on diagonals but there is no 'wrapping 
around' from row to row : the final element of each row is 
not (necessarily) equal to the first element of the subsequent 
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row. In order produce a matrix that is a BCCB matrix it 
is useful to envisage a 'virtual CCD' that is twice as big in 
both dimensions as the actual CCD, which for convenience is 
assumed to be square. This 'virtual CCD' has the property 
that it has a periodic point spread function in both direc- 
tions : the CCD is shaped like a torus. The actual CCD then 
occupies one quarter of this virtual CCD and the other three 
quarters are 'empty sky'. 

For this virtual CCD the matrix constructed by (16) 
is a {2Mf X {2Mf BCCB matrix. The first quadrant of 
this matrix, (M)^ x (M)^ in size, is identical to the matrix 
for the actual CCD. Unless explicitly stated otherwise the 
matrices for this torus-shaped virtual CCD are the ones that 
the algorithm works on. In what follows the indices 'wrap 
around' which is to say that the values of the indices are to 
be evaluated modulo the matrix dimension. 

First consider a matrix that is fully circulant rather 
than block circulant with circulant blocks. Such a matrix A 
satisfies : 



A. 



i-\-Ti j-\-7i — -^ij 



Vn, (j + n)mod(MA), 
(j + n)mod(MA) 



(Al) 



in which the matrix A has dimensions Ma x Ma- An ordi- 
nary matrix multiplication of two circulant matrices A and 
B satisfies : 



Cij — AikB, 



kj 



k 

= Ai+n iBi j+n with l = k + n 



(A2) 



Ci- 



■nj+n 



Vn 



Thus C is circulant if A and B are circulant. It is also 
straightforward to demonstrate that if A and C are circu- 
lant, B must be circulant as well : 



■ ^ ] Aj^nlBl jj^n — AilBlj 
I 

: Aik (Bfc+„ j+„ — Bkj) k = 1 — n 



(A3) 



Since A is an arbitrary circulant matrix the final equality 
can only lead to (and thus consistency) if B is indeed 
circulant. 

Since the identity matrix is itself a circulant matrix, a 
direct consequence is that the inverse of any circulant ma- 
trix A~^ is also circulant, which is trivially demonstrated by 
substituting I for C. 

If the matrix C is circulant then it is fully determined 
by its first row, as can be demonstrated by taking n = 1 — i 
in (A2), for which the following holds : 



■Y^AuB 



(A4) 



In-l-l-i 





block M, 



V 



Figure 9. The first row of blocks of a BCCB matrix ar stacked 
above each other. The first step in multiplying two BCCB matri- 
ces is to do a FT on these stacks for each matrix in the direction 
indicated by the arrow. Since the blocks individually are circu- 
lant one only needs to do this for the front face of this cube by 
performing an FT on each row and multiplying in the Fourier 
domain. 



by the following integral : 

C{x) = j dxA{x)B{ X — X ] 

which in discretized form is 
Ci = ^ viAiBi+i-i 



(A5) 



(A6) 



A one-dimensional convolution of two functions is defined 



where the vi are integration weights. The vi can always be 
arranged to be unity and so it is clear that (A6) and (A4) are 
identical summations. This means that the rrmltiplication of 
two circulant matrices can bo regarded as a convolution. 
The product of two circulant matrices can therefore be de- 
termined by multiplying the Fourier Transform (FT) of the 
first row of each of the two matrices, and then taking the 
inverse FT of this product. This yields only the first row of 
that product but since it is known to be a circulant matrix 
the other rows are then trivially found by shifting and wrap- 
ping around. Similarly the inverse of a circulant matrix can 
be found by dividing the FT of the first row of the identity 
matrix by the FT of first row of the circulant matrix to be 
inverted, and taking the inverse FT. 

If one has a BCCB matrix with dimensions x M\ in 
which the blocks have dimensions Ma x Ma , it is fully deter- 
mined by its first row of blocks. In equations (A2)-(A6) it is 
nowhere used that the individual matrix elements must be 
scalars. Thus two block circulant matrices can be multiplied 
by a one- dimensional convolution of the first row of blocks 
of each matrix, i.e. doing the same operation of convolution 
on every element of each block. One can visualize this as in 
diagram 9 by stacking the blocks and doing a one dimen- 
sional FT along columns. Since the blocks are circulant it is 
not necessary to do this operation for all columns. One does 
Ma one- dimensional FTs along the first rows of the blocks 
for each matrix and then multiplies row by row these first 
rows. 

The first step is therefore to perform the Ma one dimen- 
sional FTs of the first row of each block, for both matrices. 
The second step is to treat each of these rows as an element 
in an array, and the two arrays corresponding to matrix A 
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and B are convolved. The two-level hierarchy of one dimen- 
sional FTs can be achieved simply by one two-dimensional 
FT, in which one direction is a horizontal one and the other 
the vertical in diagram 9. The two-dimensional FT is thus 
applied to a matrix of size Ma x Ma instead of a one dimen- 
sional FT on the first row of the full matrix. For a multi- 
plication of two matrices the two-dimensional tableaux are 
multiplied element by element in the Fourier domain and 
the result is inverse FT'd. The result is the first row of the 
BCCB product matrix, where the other rows are obtained 
by shifting and wrapping around. The inverse of a matrix in 
the Fourier domain is a simple division and therefore carried 
out analogously to the multiplication of matrices. 

By using FFT algorithms the inversion of block circu- 
lant matrices with circulant blocks with size M x M is thus 
carried out in ^^(M^logM) operations. 

One more step is necessary in order to be able to apply 
this method of inverting matrices to the problem at hand. 
Because of the constraint (15) the block circulant matrix 
with circulant blocks is augmented with one row and col- 
umn. All elements of this row and column are equal to unity 
except for the corner element which is 0. Thus the matrix 
to be inverted is A' : 



A': 



(A7) 



a' o) 

where A is a BCCB matrix and 1 is a column vector of 
which all M\ elements are equal to 1. The inverse of this 
partitioned matrix can be written as (cf. Press et al., 1992) : 



A'-' = 



P Q 



(A8) 



in which P has the same dimensions as ^, Q is a column 
vector, and s is a scalar : 

P = A-'--A-' 
s 

Q = -A-'-l 
s 



l l^- A-' 



(A9) 



s 



= (l^-A-'-l) 



Since ^4 is a block circulant matrix with circulant blocks the 

column vector resulting from the product A~^ ■ 1 is itself a 
column vector al where a is some number that depends on 
A. Using this (A9) can be simplified further : 



■A- 



s = Mia 



(AlO) 



The matrix formed by 1 ■ l''" is clearly a circulant matrix 
and therefore also block circulant with circulant blocks. The 
matrix multiplications to evaluate P can therefore be done 
with FTs as described above. 

The first step in this process is sectioning the first rows 
and rearranging. The Ma x Ma matrix formed from section- 
ing and rearranging the first row of / has only one non-zero 
element which is the first element of the first row (equal to 
1). The FT of this matrix is a matrix of which every element 
is equal to 1. The matrix formed from sectioning and rear- 
ranging the first row of -Ari ■ 1^ is an Ma x Ma matrix of 



which every element is equal to and so the first element 

of the first row of its FT is equal to 1 and all other elements 
are 0. Subtracting one from the other in the Fourier domain 
produces a matrix with as the first clement of the first row 
a and all other elements equal to 1. Evaluating the FT of 
P and doing the inverse FT is then trivial. 

All the vectors b of equation (13) collected for all pixels 
also form an M'a x A4a BCCB matrix B with one extra bot- 
tom row of I's arising from the constraint (15). The product 
of this with the A'~^ is therefore also most easily carried out 
in the Fourier domain using P. Adding the final element due 
to the addition of the element Q ■ 1"^ is a simple addition of 
1/Ml to every element of P ■ i3. Since FTs are linear this 
addition can also be done in the Fourier domain, by setting 
the first element of the first row oi P ■ B equal to 1 before 
performing the inverse FT. 
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